set.seed(1)
# Data manipulation & cleaning
library(tidyverse)
library(reshape2)
library(nlme)        # gapply
library(lubridate)   # date cleaning

# Viz
library(gridExtra)
library(ggmap)
library(maps)
library(mapdata)

# Modeling & stats
library(mgcv)       # gam
library(R0)
library(Kendall)    # Mann-Kendall

# Time series & forecasting
library(forecast)
library(zoo)
library(astsa)
library(fpp2)

TODO

RT ANALYSIS

#read in full DSHS data with county and TSA information
# TOOD: rename
full.dshs <- read.csv("combined-datasets/county.csv")
tsa.dshs <- read.csv("combined-datasets/tsa.csv")
phr.df <- read.csv("combined-datasets/phr.csv")
rt.data.clean<-function(covid.data) {
  
  # select relevant metadata & DSHS population estimates
  rt.cols <- c("County", "Date", "TSA_Name", "TSA", "PHR", "PHR_Name",
               "Metro_Area", "Cases_Daily", "Population_DSHS")
  
  covid.data <- covid.data[, names(covid.data) %in% rt.cols]
             
  # set correct types
  covid.data$Date <- as.Date(covid.data$Date)
  covid.data$Cases_Daily <- as.numeric(covid.data$Cases_Daily)
  
  # TODO: decide on how to handle drops in cumulative counts -> negative daily cases
  covid.data <- covid.data %>% 
    mutate(Cases_Daily = ifelse(Cases_Daily < 0, 0, Cases_Daily))
  
  return(covid.data)
}
#separate county, metro, TSA and state data into separate dataframes
rt.data.organize <- function(mycounty, mytsa, myphr){
  
  ### COUNTY ###
  cleaned.county <- rt.data.clean(mycounty)
  county.df <- cleaned.county%>%dplyr::select(-TSA_Name, -Metro_Area, -TSA)
  
  ### TSA ###
  TSA.df <- rt.data.clean(mytsa)

  ### PHR ###
  PHR.df <- rt.data.clean(myphr)
  ### METRO ###
  # calculate daily cases and population at metro level, drop repeated rows
  metro.temp <- cleaned.county %>%
    group_by(Metro_Area,Date) %>%
    mutate(metro_cases_daily=sum(Cases_Daily, na.rm=TRUE)) %>%
    mutate(metro_pop_DSHS=sum(Population_DSHS, na.rm=TRUE)) %>%
    dplyr::select(Date, Metro_Area, metro_cases_daily, metro_pop_DSHS) %>%
    distinct()
  
  
  metro.df <- data.frame(Date = metro.temp$Date,
                         Metro_Area = metro.temp$Metro_Area,
                         Cases_Daily = metro.temp$metro_cases_daily,
                         Population_DSHS = metro.temp$metro_pop_DSHS)
  
  
  ### STATE ###
  # calculate daily cases and population at state level, drop repeated rows
  state.temp <- TSA.df %>%
    group_by(Date) %>%
    mutate(state_cases_daily=sum(Cases_Daily, na.rm=TRUE)) %>%
    mutate(state_pop_DSHS=sum(Population_DSHS, na.rm=TRUE)) %>%
    dplyr::select(Date, state_cases_daily, state_pop_DSHS) %>%
    distinct()
  
  
  state.df <- data.frame(Date = state.temp$Date, 
                         Cases_Daily = state.temp$state_cases_daily,
                         Population_DSHS = state.temp$state_pop_DSHS)

  ### OUTPUT ###
  rt.df.out <- list(county=county.df, TSA=TSA.df, state=state.df, metro=metro.df, phr = PHR.df)
  return(rt.df.out)
}
rt.df.extraction <- function(Rt.estimate.output) {

  # extract r0 estimate values into dataframe
  rt.df <- setNames(stack(Rt.estimate.output$estimates$TD$R)[2:1], c('Date', 'Rt'))
  rt.df$Date <- as.Date(rt.df$Date)

  # get 95% CI
  CI.lower.list <- Rt.estimate.output$estimates$TD$conf.int$lower
  CI.upper.list <- Rt.estimate.output$estimates$TD$conf.int$upper

  #use unlist function to format as vector
  CI.lower <- unlist(CI.lower.list, recursive = TRUE, use.names = TRUE)
  CI.upper <- unlist(CI.upper.list, recursive = TRUE, use.names = TRUE)

  rt.df$lower <- CI.lower
  rt.df$upper <- CI.upper
  
  rt.df <- rt.df %>%
    mutate(lower = replace(lower, Rt == 0, NA)) %>%
    mutate(upper = replace(upper, Rt == 0, NA)) %>%
    mutate(Rt = replace(Rt, Rt == 0, NA))
  
  return(rt.df)
}
#Function to generate Rt estimates, need to read in data frame and population size
covid.rt <- function(mydata) {
  
  ### DECLARE VALS ###
  
  #set generation time
  #Tapiwa, Ganyani "Esimating the gen interval for Covid-19":

  # LOGNORMAL OPS
  # gen.time<-generation.time("lognormal", c(4.0, 2.9))
  # gen.time<-generation.time("lognormal", c(4.7,2.9)) #Nishiura

  # GAMMA OPS
  # gen.time<-generation.time("gamma", c(5.2, 1.72)) #Singapore
  # gen.time<-generation.time("gamma", c(3.95, 1.51)) #Tianjin
  gen.time <- generation.time("gamma", c(3.96, 4.75))
  print(as.character(mydata[1,2]))
  
  # TODO: consider removing this and handling na cases directly
  #change na values to 0
  mydata <- mydata %>% replace(is.na(.),0)
  sum.daily.cases <- sum(mydata$Cases_Daily)
  pop.DSHS <- mydata$Population_DSHS[1]
  
  #Get 7 day moving average of daily cases
  mydata$MA_7day<-rollmean(mydata$Cases_Daily, k=7, na.pad=TRUE, align='right')
  #change na values to 0
  mydata<-mydata%>%replace(is.na(.),0)

  #create a vector of new cases 7 day moving average
  mydata.new<-pull(mydata, MA_7day)
  #mydata.new <- pull(mydata, Cases_Daily)
  
  # get dates as vectors
  date.vector <- pull(mydata, Date)
  
  #create a named numerical vector using the date.vector as names of new cases
  #Note: this is needed to run R0 package function estimate.R()
  names(mydata.new) <- c(date.vector)
  
  #Find min.date, max.date and row number of max.date
  min.date <- min(mydata$Date)
  max.date <- max(mydata$Date)

  #get row number of March 15 and first nonzero entry
  #NOTE: for 7 day moving average start March 15, for daily start March 9
  #find max row between the two (this will be beginning of rt data used)
  march15.row <- which(mydata$Date=="2020-03-15")
  first.nonzero <- min(which(mydata$Cases_Daily>0))
  last.nonzero <- max(which(mydata$Cases_Daily>0))
  minrow <- max(march15.row, first.nonzero)

  ### R0 ESTIMATION ###
  # TODO: establish better criteria for # of required daily cases
  # mean daily cases; average last x number of days; % of region pop etc.
  if(!is.na(minrow) & sum.daily.cases > 100) {
    
    # declare limits for rt estimation (first nonzero date & last nonzero date)
    i <- minrow
    j <- as.integer(min(last.nonzero, max.date))

    #reduce the vector to be rows from min date (March 9 or first nonzero case) to current date
    mydata.newest <- mydata.new[i:j]
    rt.DSHS <- estimate.R(mydata.newest, 
                        gen.time, 
                        begin=as.integer(1),
                        end=length(mydata.newest),
                        methods=c("TD"), 
                        pop.size=pop.DSHS,
                        nsim=1000)
  
    rt.DSHS.df <- rt.df.extraction(rt.DSHS)
    
  } else {
    # error catch small regions & na values for minrow
    rt.DSHS.df <- data.frame(Date = as.Date(mydata$Date),
                             Rt = rep(NA, length(mydata$Date)),
                             lower = rep(NA, length(mydata$Date)),
                             upper = rep(NA, length(mydata$Date)))
  }
  
  return(rt.DSHS.df)
}

County

#apply the covid.rt function to the county.datily data by county 
county.rt.output <- nlme::gapply(county.daily, FUN = covid.rt, groups = county.daily$County)
## [1] "Anderson"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Andrews"
## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.142857142857143, `2020-04-05` =
## 0.857142857142857, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.142857142857143, `2020-04-05` =
## 0.857142857142857, : Using initial incidence as initial number of cases.
## [1] "Angelina"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Aransas"
## [1] "Archer"
## [1] "Armstrong"
## [1] "Atascosa"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Austin"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Bailey"
## Warning in est.R0.TD(epid = c(`2020-05-05` = 0.142857142857143, `2020-05-06` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-05-05` = 0.142857142857143, `2020-05-06` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Bandera"
## [1] "Bastrop"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Baylor"
## [1] "Bee"
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Bell"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Bexar"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Blanco"
## [1] "Borden"
## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf
## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf
## [1] "Bosque"
## [1] "Bowie"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Brazoria"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Brazos"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Brewster"
## Warning in est.R0.TD(epid = c(`2020-05-01` = 0.142857142857143, `2020-05-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-05-01` = 0.142857142857143, `2020-05-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Briscoe"
## [1] "Brooks"
## [1] "Brown"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Burleson"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Burnet"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Caldwell"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Calhoun"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Callahan"
## [1] "Cameron"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Camp"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Carson"
## [1] "Cass"
## [1] "Castro"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Chambers"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Cherokee"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Childress"
## [1] "Clay"
## [1] "Cochran"
## [1] "Coke"
## [1] "Coleman"
## [1] "Collin"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Using initial incidence as initial number of cases.
## [1] "Collingsworth"
## [1] "Colorado"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.285714285714286, `2020-04-02` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.285714285714286, `2020-04-02` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Comal"
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.428571428571429, `2020-03-23` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.428571428571429, `2020-03-23` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Comanche"
## [1] "Concho"
## [1] "Cooke"
## Warning in est.R0.TD(epid = c(`2020-04-10` = 0.142857142857143, `2020-04-11` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-10` = 0.142857142857143, `2020-04-11` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Coryell"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Cottle"
## [1] "Crane"
## [1] "Crockett"
## [1] "Crosby"
## [1] "Culberson"
## [1] "Dallam"
## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Dallas"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 1.14285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 1.14285714285714, : Using initial incidence as initial number of cases.
## [1] "Dawson"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Deaf Smith"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Delta"
## [1] "Denton"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "DeWitt"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.142857142857143, `2020-03-20` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.142857142857143, `2020-03-20` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Dickens"
## [1] "Dimmit"
## [1] "Donley"
## [1] "Duval"
## Warning in est.R0.TD(epid = c(`2020-04-14` = 0.142857142857143, `2020-04-15` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-14` = 0.142857142857143, `2020-04-15` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Eastland"
## [1] "Ector"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Edwards"
## [1] "El Paso"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Ellis"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Erath"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Falls"
## [1] "Fannin"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Fayette"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Fisher"
## [1] "Floyd"
## [1] "Foard"
## [1] "Fort Bend"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Franklin"
## [1] "Freestone"
## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Frio"
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Gaines"
## [1] "Galveston"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Garza"
## [1] "Gillespie"
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Glasscock"
## [1] "Goliad"
## [1] "Gonzales"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Gray"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Grayson"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Gregg"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Grimes"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Guadalupe"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.571428571428571, `2020-03-26` =
## 1.14285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.571428571428571, `2020-03-26` =
## 1.14285714285714, : Using initial incidence as initial number of cases.
## [1] "Hale"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hall"
## [1] "Hamilton"
## [1] "Hansford"
## [1] "Hardeman"
## [1] "Hardin"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Harris"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Using initial incidence as initial number of cases.
## [1] "Harrison"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hartley"
## [1] "Haskell"
## [1] "Hays"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hemphill"
## [1] "Henderson"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hidalgo"
## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.285714285714286, `2020-03-24` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.285714285714286, `2020-03-24` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Hill"
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hockley"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Hood"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Hopkins"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Houston"
## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Howard"
## [1] "Hudspeth"
## [1] "Hunt"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hutchinson"
## [1] "Irion"
## [1] "Jack"
## [1] "Jackson"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Jasper"
## [1] "Jeff Davis"
## [1] "Jefferson"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 1, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 1, : Using initial incidence as initial number of cases.
## [1] "Jim Hogg"
## [1] "Jim Wells"
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Johnson"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Jones"
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Karnes"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Kaufman"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Kendall"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Kenedy"
## [1] "Kent"
## [1] "Kerr"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Kimble"
## [1] "King"
## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf
## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf
## [1] "Kinney"
## [1] "Kleberg"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Knox"
## [1] "La Salle"
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Lamar"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Lamb"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Lampasas"
## [1] "Lavaca"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Lee"
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Leon"
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Liberty"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Limestone"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Lipscomb"
## [1] "Live Oak"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Llano"
## [1] "Loving"
## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf
## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf
## [1] "Lubbock"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Lynn"
## [1] "Madison"
## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Marion"
## [1] "Martin"
## [1] "Mason"
## [1] "Matagorda"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Maverick"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "McCulloch"
## [1] "McLennan"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Using initial incidence as initial number of cases.
## [1] "McMullen"
## [1] "Medina"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Menard"
## [1] "Midland"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Milam"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Mills"
## [1] "Mitchell"
## [1] "Montague"
## [1] "Montgomery"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Moore"
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Morris"
## [1] "Motley"
## [1] "Nacogdoches"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Navarro"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Newton"
## [1] "Nolan"
## Warning in est.R0.TD(epid = c(`2020-04-24` = 0.142857142857143, `2020-04-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-24` = 0.142857142857143, `2020-04-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Nueces"
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Ochiltree"
## [1] "Oldham"
## [1] "Orange"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Palo Pinto"
## [1] "Panola"
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.428571428571429, `2020-04-03` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.428571428571429, `2020-04-03` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "Parker"
## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.142857142857143, `2020-03-24` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.142857142857143, `2020-03-24` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Parmer"
## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Pecos"
## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Polk"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Potter"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.285714285714286, `2020-03-22` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.285714285714286, `2020-03-22` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Presidio"
## [1] "Rains"
## [1] "Randall"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.428571428571429, `2020-03-29` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.428571428571429, `2020-03-29` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Reagan"
## [1] "Real"
## [1] "Red River"
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Reeves"
## [1] "Refugio"
## [1] "Roberts"
## [1] "Robertson"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Rockwall"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Runnels"
## [1] "Rusk"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Sabine"
## [1] "San Augustine"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "San Jacinto"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "San Patricio"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "San Saba"
## [1] "Schleicher"
## [1] "Scurry"
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Shackelford"
## [1] "Shelby"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Sherman"
## [1] "Smith"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Somervell"
## [1] "Starr"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Stephens"
## [1] "Sterling"
## [1] "Stonewall"
## [1] "Sutton"
## [1] "Swisher"
## [1] "Tarrant"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Taylor"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Terrell"
## [1] "Terry"
## [1] "Throckmorton"
## [1] "Titus"
## Warning in est.R0.TD(epid = c(`2020-04-03` = 0.142857142857143, `2020-04-04` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-03` = 0.142857142857143, `2020-04-04` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Tom Green"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Travis"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Trinity"
## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.285714285714286, `2020-04-05` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.285714285714286, `2020-04-05` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Tyler"
## [1] "Upshur"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Upton"
## [1] "Uvalde"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Val Verde"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Van Zandt"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Victoria"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.428571428571429, `2020-03-27` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.428571428571429, `2020-03-27` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Walker"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Waller"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.428571428571429, `2020-03-30` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.428571428571429, `2020-03-30` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Ward"
## [1] "Washington"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.714285714285714, `2020-03-29` =
## 0.714285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.714285714285714, `2020-03-29` =
## 0.714285714285714, : Using initial incidence as initial number of cases.
## [1] "Webb"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Wharton"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Wheeler"
## [1] "Wichita"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Wilbarger"
## [1] "Willacy"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Williamson"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 1.14285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 1.14285714285714, : Using initial incidence as initial number of cases.
## [1] "Wilson"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Winkler"
## [1] "Wise"
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Wood"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Yoakum"
## [1] "Young"
## [1] "Zapata"
## Warning in est.R0.TD(epid = c(`2020-04-07` = 0.142857142857143, `2020-04-08` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-07` = 0.142857142857143, `2020-04-08` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Zavala"
#convert list of dataframes to single dataframe with county as id col
county.rt.df <- data.table::rbindlist(county.rt.output, idcol = TRUE)
colnames(county.rt.df)[1] = 'County'

Metro

#apply the covid.rt function to the metro.daily data by metro region
metro.rt.output <- nlme::gapply(metro.daily, FUN=covid.rt, groups=metro.daily$Metro_Area)
## [1] "Metro"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 7, `2020-03-16` = 6, `2020-03-17` =
## 6.57142857142857, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 7, `2020-03-16` = 6, `2020-03-17` =
## 6.57142857142857, : Using initial incidence as initial number of cases.
## [1] "Non-Metro"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
#convert list of data frames to single data frame with metro as id col
metro.rt.df <- data.table::rbindlist(metro.rt.output, idcol=TRUE)
colnames(metro.rt.df)[1] = 'Metro_Area'

TSA

# plot with shaded confidence intervals
rt.plot<-function(rt.data, caption){
  library(ggplot2)
    caption<-rt.data$TSA[1]
  plot<-ggplot(rt.data, aes(Date, Rt))+
    geom_ribbon(aes(ymin=lower, ymax=upper), fill="light grey")+
    geom_point(color="blue", size=0.2)+
    labs(title=caption)
  plot
}
#apply the covid.rt function to the TSA.daily data by TSA region
TSA.rt.output <- nlme::gapply(TSA.daily, FUN=covid.rt, groups=TSA.daily$TSA)
## [1] "A"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.571428571428571, `2020-03-22` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.571428571428571, `2020-03-22` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "B"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.714285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.714285714285714, : Using initial incidence as initial number of cases.
## [1] "C"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "D"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "E"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Using initial incidence as initial number of cases.
## [1] "F"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "G"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "H"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "I"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "J"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "K"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "L"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "M"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Using initial incidence as initial number of cases.
## [1] "N"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "O"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "P"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Q"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.57142857142857, `2020-03-16` =
## 1.71428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.57142857142857, `2020-03-16` =
## 1.71428571428571, : Using initial incidence as initial number of cases.
## [1] "R"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "S"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "T"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "U"
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "V"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
#convert list of data frames to single data frame with TSA as variable
TSA.rt.df <- data.table::rbindlist(TSA.rt.output, idcol=TRUE)

# add TSA_Name & fix dates
colnames(TSA.rt.df)[1] <- 'TSA'
TSA.rt.df$Date <- as.Date(TSA.rt.df$Date)
TSA.rt.df <- merge(TSA.rt.df, TSA.daily[, c('TSA', 'TSA_Name', 'Date')], by = c('TSA', 'Date'))

# combine TSA 
TSA.rt.df$TSA <- paste0(TSA.rt.df$TSA, ' - ', TSA.rt.df$TSA_Name)
TSA.rt.df$TSA_Name <- NULL

#plot Rt by TSA
nlme::gapply(TSA.rt.df, FUN=rt.plot, groups=TSA.rt.df$TSA)
## $`A - Amarillo`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`B - Lubbock`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`C - Wichita Falls`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`D - Abilene`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`E - Dallas/Ft. Worth`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`F - Paris`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`G - Longview/Tyler`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`H - Lufkin`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`I - El Paso`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`J - Midland/Odessa`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`K - San Angelo`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`L - Belton/Killeen`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`M - Waco`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`N - Bryan/College Station`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`O - Austin`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`P - San Antonio`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`Q - Houston`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`R - Galveston`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`S - Victoria`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`T - Laredo`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`U - Corpus Christi`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`V - Lower Rio Grande Valley`
## Warning: Removed 1 rows containing missing values (geom_point).

PHR

phr.rt.output <- nlme::gapply(phr.daily, FUN=covid.rt, groups=phr.daily$PHR)
## [1] "1"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 1.28571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 1.28571428571429, : Using initial incidence as initial number of cases.
## [1] "11"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "2/3"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Using initial incidence as initial number of cases.
## [1] "4/5N"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "6/5S"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 3, `2020-03-16` =
## 2.14285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 3, `2020-03-16` =
## 2.14285714285714, : Using initial incidence as initial number of cases.
## [1] "7"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "8"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "9/10"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
#convert list of data frames to single data frame with phr as id col
phr.rt.df <- data.table::rbindlist(phr.rt.output, idcol=TRUE)

# add PHR_Name & fix dates
colnames(phr.rt.df)[1] <- 'PHR'
phr.rt.df$Date <- as.Date(phr.rt.df$Date)
phr.rt.df <- merge(phr.rt.df, phr.daily[, c('PHR', 'PHR_Name', 'Date')], by = c('PHR', 'Date'))

# combine phr 
phr.rt.df$PHR <- paste0(phr.rt.df$PHR, ' - ', phr.rt.df$PHR_Name)
phr.rt.df$PHR_Name <- NULL

State

state.rt.df <- covid.rt(state.daily)
## [1] "0"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 7.28571428571429, `2020-03-16` =
## 6.28571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 7.28571428571429, `2020-03-16` =
## 6.28571428571429, : Using initial incidence as initial number of cases.

export

write.csv(county.rt.df,"statistical-output/rt/county_rt.csv", row.names = FALSE)
write.csv(metro.rt.df, "statistical-output/rt/metro_rt.csv", row.names = FALSE)
write.csv(TSA.rt.df, "statistical-output/rt/tsa_rt.csv", row.names = FALSE)
write.csv(phr.rt.df, "statistical-output/rt/phr_rt.csv", row.names = FALSE)
write.csv(state.rt.df, "statistical-output/rt/state_rt.csv", row.names = FALSE)

grouping

colnames(county.rt.df)[1] <- 'Level'
colnames(metro.rt.df)[1] <- 'Level'
colnames(TSA.rt.df)[1] <- 'Level'
colnames(phr.rt.df)[1] <- 'Level'
state.rt.df$Level <- 'Texas'

county.rt.df$Level_Type = 'County'
metro.rt.df$Level_Type = 'Metro'
TSA.rt.df$Level_Type = 'TSA'
phr.rt.df$Level_Type = 'PHR'
state.rt.df$Level_Type = 'State'

combined.rt.df <- rbind(subset(county.rt.df, Date != max(Date)),
                        subset(metro.rt.df, Date != max(Date)),
                        subset(TSA.rt.df, Date != max(Date)),
                        subset(state.rt.df, Date != max(Date)),
                        subset(phr.rt.df, Date != max(Date)))

write.csv(combined.rt.df, 'statistical-output/rt/stacked_rt.csv', row.names = FALSE)
write.csv(combined.rt.df, 'tableau/stacked_rt.csv', row.names = FALSE)

TIME SERIES

#ts.rt.data.clean function will clean data, dropping unwanted variables
ts.data.clean<-function(covid.data) {

  ts_cols <- c("County", "Date", "TSA", "TSA_Name", "PHR", "PHR_Name",
               "Metro_Area", "Cases_Daily", "Cases_Cumulative",
               "Population_DSHS", "Hospitalizations_Total")
  
  covid.data <- covid.data[, names(covid.data) %in% ts_cols]
  
  # set correct types
  covid.data$Date <- as.Date(covid.data$Date)
  covid.data$Cases_Daily <- as.numeric(covid.data$Cases_Daily)
  
  #change any Cases_Daily below zero to zero
  covid.data <- covid.data %>% mutate(Cases_Daily = ifelse(Cases_Daily < 0, 0, Cases_Daily))
  return(covid.data)
}
#separate county, metro, TSA and State data into separate Data frames
ts.data.organize<-function(mycounty, mytsa, myphr){

  ### COUNTY ###
  # clean county vals and restrict to first date of collection
  cleaned.county <- ts.data.clean(mycounty)
  cleaned.county <- subset(cleaned.county, !is.na(Date) & Date >= as.Date('2020-03-04'))

  # Set column names
  county.df <- cleaned.county%>%dplyr::select(-TSA_Name, -Metro_Area)
  county.df$Level <- 'County'
  colnames(county.df)[2] = 'Level_Name'
  
  #change hospitalizations to numeric
  mytsa$Hospitalizations_Total<-as.numeric(mytsa$Hospitalizations_Total)
  
  
  # PHR 
  PHR.df <- ts.data.clean(myphr)
  PHR.df$Level = 'PHR'
  colnames(PHR.df)[3] = 'Level_Name' 
  
  ### METRO ###
  metro.temp<-cleaned.county %>%
    group_by(Metro_Area, Date) %>%
    mutate(metro_Cases_Daily = sum(Cases_Daily, na.rm = TRUE)) %>%
    mutate(metro_Cases_Cumulative = sum(Cases_Cumulative, na.rm=TRUE)) %>%
    mutate(metro_pop_DSHS = sum(Population_DSHS, na.rm=TRUE)) %>%
    dplyr::select(Date, Metro_Area, metro_Cases_Daily, metro_Cases_Cumulative, metro_pop_DSHS) %>%
    distinct()
    
  metro.df <- data.frame(Date = metro.temp$Date,
                         Level = 'metro',
                         Level_Name = metro.temp$Metro_Area,
                         Cases_Daily = metro.temp$metro_Cases_Daily,
                         Cases_Cumulative = metro.temp$metro_Cases_Cumulative,
                         Population_DSHS = metro.temp$metro_pop_DSHS)
  
  # drop NA dates
  metro.df <- subset(metro.df, !is.na(Date) & Date>=as.Date('2020-03-04'))
  
  ### TSA ###
  TSA.df <- ts.data.clean(mytsa)
  TSA.df <- subset(TSA.df, !is.na(Date) & Date >= as.Date('2020-03-04'))
  
  TSA.df$Level <- 'TSA'
  colnames(TSA.df)[2] <- 'Level_Name'
  
  
  ### STATE ###
  state.temp<-TSA.df %>%
    group_by(Date) %>%
    mutate(state_Cases_Daily = sum(Cases_Daily, na.rm = TRUE)) %>%
    mutate(state_Cases_Cumulative = sum(Cases_Cumulative, na.rm = TRUE)) %>%
    mutate(state_pop_DSHS = sum(Population_DSHS, na.rm = TRUE)) %>%
    mutate(state_hosp = sum(Hospitalizations_Total, na.rm=TRUE))%>%
    dplyr::select(Date, state_Cases_Daily, state_Cases_Cumulative, state_pop_DSHS, state_hosp) %>%
    distinct()
    
  state.df <- data.frame(Date=state.temp$Date,
                         Level = 'State',
                         Level_Name = 'Texas',
                         Cases_Daily=state.temp$state_Cases_Daily,
                         Cases_Cumulative=state.temp$state_Cases_Cumulative, 
                         Population_DSHS=state.temp$state_pop_DSHS,
                         Hospitalizations_Total=state.temp$state_hosp)
  
  state.df <- subset(state.df, Date>=as.Date('2020-03-04'))
  
  
  ### OUTPUT ###
  ts.df.out <- list(county=county.df, TSA=TSA.df, state=state.df, metro=metro.df, phr = PHR.df)
  return(ts.df.out)
}
# Compute forecast (UPDATE PREDICTION PERIOD [days] AS NEEDED)
covid.arima.forecast<-function(mydata, prediction.period = 10, mindate)
{
  maxdate <- max(mydata$Date)
  # mindate <- as.Date('2020-05-01')
  pred_start_label = format(mindate, format = '%m_%d')
  
  mydata = subset(mydata, Date >= mindate)
  model.length <- as.numeric(length(mydata$Date) + prediction.period)

  if(max(mydata$Cases_Daily>=100, na.rm = TRUE))
  {
    # arima requires cases to be a timeseries vector
    # convert daily cases to time series
    my.timeseries<-ts(mydata$Cases_Daily)
    print(mydata[1,2])
    #load package(pracma)
    library(pracma)
    my.timeseries<-movavg(my.timeseries,7,"s")
    #d=0 restricts first differencing to 0 so that daily cases aren't differenced
    
    arima.fit <- forecast::auto.arima(my.timeseries)
  
    # obtain diagnostic plots for ideal arima (p,d,q) selection 
    acf <- ggAcf(my.timeseries, lag.max=30)
    pacf <- ggPacf(my.timeseries, lag.max=30)
    ts.diagnostics <- grid.arrange(acf, pacf, nrow=2)
    ggsave(paste0('statistical-output/time-series/diagnostics/',
                  mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'),
           plot = ts.diagnostics)

    
    
    # save parameters from arima autofit
    p<- arima.fit$arma[1]          # autoregressive order 
    q<- arima.fit$arma[2]          # moving average order 
    d<-arima.fit$arma[6]           # differencing order from model
  
    # 10 day forecast, CI for lower and upper has confidence level 70% set by level =c(70,70)
    arima.forecast <- forecast::forecast(arima.fit, h = prediction.period, level=c(70,70))

    #return a dataframe of the arima model(Daily cases by date)
    arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
                            Cases_Raw = c(mydata$Cases_Daily, rep(NA, times = prediction.period)),
                            Cases_Daily = c(my.timeseries, arima.forecast[['mean']]),
                            CI_Lower = c(rep(NA, times = length(my.timeseries)),
                                         arima.forecast[['lower']][, 2]),
                            CI_Upper = c(rep(NA, times = length(my.timeseries)),
                                         arima.forecast[['upper']][, 2]))
                            # Order_AutoReg = c(rep(p, times = model.length)),
                            # Order_Moving_Avg = c(rep(q, times = model.length)),
                            # Differencing = c(rep(d, times = model.length)))
    
      # save prediction plot for preliminary review
      arima.plot <- ggplot(arima.out, aes(x=Date, y = Cases_Daily))+
                    geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), fill = "red", alpha = 0.5, size = 0.1)+
                    geom_line(color = "black", size = 1)+
                    labs(y = 'Daily Cases (7-Day Moving Average)', x = 'Date',
                         title = mydata$Level_Name[1]) +
                    scale_x_date(limits = c(mindate, maxdate + prediction.period),
                                 date_labels = '%b-%d', date_breaks = '1 week') + 
                    ggpubr::theme_pubr() +
                    theme(axis.text.x = element_text(angle = -90))
      
      ggsave(plot = arima.plot, paste0('statistical-output/time-series/plots/',
                                       mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'))    
    
    } else {
    # insufficient data catch: return NA values for predictions 
    arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
                            Cases_Raw = c(mydata$Cases_Daily, rep(NA, times = prediction.period)),
                            Cases_Daily = rep(NA, times = model.length),
                            CI_Lower = rep(NA, times = model.length),
                            CI_Upper =  rep(NA, times = model.length))
                            # Order_AutoReg = c(rep(NA, times = model.length)),
                            # Order_Moving_Avg = c(rep(NA, times = model.length)),
                            # Differencing = c(rep(NA, times= model.length)))
    }
  
  #replace CI lower limit with 0 if negative
  arima.out$CI_Lower <- ifelse(arima.out$CI_Lower>=0,arima.out$CI_Lower, 0)
  
  
  return(arima.out)
}

County

# apply arima across all counties
county.arima.output.1 <- nlme::gapply(county.Daily,
                                    FUN = covid.arima.forecast,
                                    groups = county.Daily$Level_Name,
                                    mindate = as.Date('2020-03-04'))
## [1] "Anderson"
## 
## Attaching package: 'pracma'
## The following object is masked from 'package:mgcv':
## 
##     magic
## The following object is masked from 'package:purrr':
## 
##     cross
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Angelina"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Bell"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Bexar"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Brazoria"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Brazos"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Cameron"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Chambers"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Collin"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Comal"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Dallas"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Denton"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Ector"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "El Paso"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Ellis"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Fort Bend"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Galveston"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Gregg"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Grimes"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Harris"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Hays"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Hidalgo"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Jefferson"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Johnson"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Jones"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Kaufman"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "La Salle"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Lubbock"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Maverick"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "McLennan"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Medina"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Midland"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Montgomery"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Moore"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Nueces"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Orange"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Potter"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Randall"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Rusk"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Scurry"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Smith"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Starr"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Tarrant"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Titus"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Tom Green"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Travis"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Val Verde"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Victoria"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Walker"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Webb"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Williamson"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Wilson"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

# bind list of dataframes to one dataframe
county.arima.df.1 <- data.table::rbindlist(county.arima.output.1, idcol = TRUE)
colnames(county.arima.df.1)[1] <- 'County'

# county.arima.output.2 <- nlme::gapply(county.Daily,
#                                     FUN = covid.arima.forecast,
#                                     groups = county.Daily$Level_Name,
#                                     mindate = as.Date('2020-06-03'))
# 
# 
# # bind list of dataframes to one dataframe
# county.arima.df.2 <- data.table::rbindlist(county.arima.output.2, idcol = TRUE)
# colnames(county.arima.df.2)[1] <- 'County'

Metro

# apply arima across both metro values
metro.arima.output.1 <- nlme::gapply(metro.Daily,
                                   FUN = covid.arima.forecast,
                                   groups = metro.Daily$Level_Name,
                                   mindate = as.Date('2020-03-04'))
## [1] "metro"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "metro"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

# bind list of dataframes to one dataframe
metro.arima.df.1 <- data.table::rbindlist(metro.arima.output.1, idcol = TRUE)
colnames(metro.arima.df.1)[1] <- 'Metro_Area'
# 
# metro.arima.output.2 <- nlme::gapply(metro.Daily,
#                                    FUN = covid.arima.forecast,
#                                    groups = metro.Daily$Level_Name,
#                                    mindate = as.Date('2020-06-03'))
# 
# # bind list of dataframes to one dataframe
# metro.arima.df.2 <- data.table::rbindlist(metro.arima.output.2, idcol = TRUE)
# colnames(metro.arima.df.2)[1] <- 'Metro_Area'

TSA

# All data
TSA.arima.output.1 <- nlme::gapply(TSA.Daily,
                                   FUN = covid.arima.forecast,
                                   groups = TSA.Daily$Level_Name,
                                   mindate = as.Date('2020-03-04'))
## [1] "A"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "B"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "D"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "E"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "F"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "G"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "H"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "I"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "J"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "K"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "L"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "M"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "N"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "O"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "P"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "Q"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "R"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "S"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "T"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "U"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "V"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

TSA.arima.df.1 <- data.table::rbindlist(TSA.arima.output.1, idcol=TRUE)

colnames(TSA.arima.df.1)[1] <- 'TSA'
TSA.arima.df.1 <- merge(TSA.arima.df.1, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = c('TSA'))

# combine TSA 
TSA.arima.df.1$TSA <- paste0(TSA.arima.df.1$TSA, ' - ', TSA.arima.df.1$TSA_Name)
TSA.arima.df.1$TSA_Name <- NULL

PHR

PHR.arima.output <- nlme::gapply(PHR.Daily,
                                 FUN = covid.arima.forecast,
                                 groups = PHR.Daily$Level_Name,
                                 mindate = as.Date('2020-03-04'))
## [1] "2/3"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "9/10"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "11"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "6/5S"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "1"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "8"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "7"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] "4/5N"
## Saving 7 x 5 in image
## Saving 7 x 5 in image

PHR.arima.df <- data.table::rbindlist(PHR.arima.output, idcol=TRUE)
colnames(PHR.arima.df)[1] <- 'Level_Name'

PHR.arima.df <- merge(PHR.arima.df, unique(PHR.Daily[, c('PHR', 'Level_Name')]), by = c('Level_Name'))

# combine PHR 
PHR.arima.df$Level_Name <- paste0(PHR.arima.df$PHR, ' - ', PHR.arima.df$Level_Name)
PHR.arima.df$PHR <- NULL
colnames(PHR.arima.df)[1] <- 'PHR'

State

texas.arima.df.1 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-03-04'))
## Saving 7 x 5 in image
## Saving 7 x 5 in image

# texas.arima.df.2 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-06-03'))
#Create plot function, forecast.data is the dataframe containing true data
# and the forecast data. days.forecast is the number of days that are forecast
# in the dataframe (here we are forecasting 10 days, so it should be 10 unless we
# change this at some point down the road)
forecast.plot<-function(forecast.data, plot.title, y.label){
    mindate<-as.POSIXct("2020-03-04")
    maxdate<-as.POSIXct(max(forecast.data$Date))
  ts.plot<-ggplot(aes(x=as.POSIXct(Date), y=Cases_Daily), data=forecast.data)+
    geom_ribbon(aes(ymin=CI_Lower, ymax=CI_Upper), fill="grey50", size=0.1)+
    geom_line(color="blue", size=1)+
    geom_line(aes(x=as.POSIXct(Date), y=CI_Lower), color="grey", size=0.1)+
    geom_line(aes(x=as.POSIXct(Date),y=CI_Upper), color="grey", size=0.1)+
    scale_x_datetime(limits = c(mindate, maxdate))+
    xlab("Date")+
    ylab(y.label)+
    #Can use the following title if we are running using nlme for data frame
    #ggtitle(paste(forecast.data$.id,": TS Daily Cases + 10 Day Forcast",sep=""))
    ggtitle(plot.title)
ts.plot
}
######### View Output Table & Graph for TSA Q ##########
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
#subset TSA Q arima data
TSA_Q.arima.df<-TSA.arima.df.1%>%subset(TSA.arima.df.1$TSA== "Q - Houston")
#apply the plot function to TSA Q arima data frame
TSA_Q.arima.df%>%kable(caption="Greater Houston ARIMA - Daily Cases")%>%kable_styling(full_width=FALSE)
Greater Houston ARIMA - Daily Cases
TSA Date Cases_Raw Cases_Daily CI_Lower CI_Upper
Q - Houston 2020-03-04 0 0.000000 NA NA
Q - Houston 2020-03-05 0 0.000000 NA NA
Q - Houston 2020-03-06 5 1.666667 NA NA
Q - Houston 2020-03-07 NA NA NA NA
Q - Houston 2020-03-08 NA NA NA NA
Q - Houston 2020-03-09 6 NA NA NA
Q - Houston 2020-03-10 0 NA NA NA
Q - Houston 2020-03-11 0 NA NA NA
Q - Houston 2020-03-12 3 NA NA NA
Q - Houston 2020-03-13 0 NA NA NA
Q - Houston 2020-03-14 NA NA NA NA
Q - Houston 2020-03-15 9 NA NA NA
Q - Houston 2020-03-16 0 NA NA NA
Q - Houston 2020-03-17 0 NA NA NA
Q - Houston 2020-03-18 0 NA NA NA
Q - Houston 2020-03-19 2 NA NA NA
Q - Houston 2020-03-20 17 NA NA NA
Q - Houston 2020-03-21 2 4.285714 NA NA
Q - Houston 2020-03-22 1 3.142857 NA NA
Q - Houston 2020-03-23 0 3.142857 NA NA
Q - Houston 2020-03-24 60 11.714286 NA NA
Q - Houston 2020-03-25 109 27.285714 NA NA
Q - Houston 2020-03-26 72 37.285714 NA NA
Q - Houston 2020-03-27 60 43.428571 NA NA
Q - Houston 2020-03-28 48 50.000000 NA NA
Q - Houston 2020-03-29 242 84.428571 NA NA
Q - Houston 2020-03-30 102 99.000000 NA NA
Q - Houston 2020-03-31 52 97.857143 NA NA
Q - Houston 2020-04-01 190 109.428571 NA NA
Q - Houston 2020-04-02 227 131.571429 NA NA
Q - Houston 2020-04-03 138 142.714286 NA NA
Q - Houston 2020-04-04 212 166.142857 NA NA
Q - Houston 2020-04-05 204 160.714286 NA NA
Q - Houston 2020-04-06 138 165.857143 NA NA
Q - Houston 2020-04-07 473 226.000000 NA NA
Q - Houston 2020-04-08 477 267.000000 NA NA
Q - Houston 2020-04-09 231 267.571429 NA NA
Q - Houston 2020-04-10 781 359.428571 NA NA
Q - Houston 2020-04-11 275 368.428571 NA NA
Q - Houston 2020-04-12 333 386.857143 NA NA
Q - Houston 2020-04-13 76 378.000000 NA NA
Q - Houston 2020-04-14 163 333.714286 NA NA
Q - Houston 2020-04-15 292 307.285714 NA NA
Q - Houston 2020-04-16 236 308.000000 NA NA
Q - Houston 2020-04-17 267 234.571429 NA NA
Q - Houston 2020-04-18 277 234.857143 NA NA
Q - Houston 2020-04-19 211 217.428571 NA NA
Q - Houston 2020-04-20 194 234.285714 NA NA
Q - Houston 2020-04-21 216 241.857143 NA NA
Q - Houston 2020-04-22 201 228.857143 NA NA
Q - Houston 2020-04-23 195 223.000000 NA NA
Q - Houston 2020-04-24 181 210.714286 NA NA
Q - Houston 2020-04-25 214 201.714286 NA NA
Q - Houston 2020-04-26 171 196.000000 NA NA
Q - Houston 2020-04-27 161 191.285714 NA NA
Q - Houston 2020-04-28 139 180.285714 NA NA
Q - Houston 2020-04-29 219 182.857143 NA NA
Q - Houston 2020-04-30 276 194.428571 NA NA
Q - Houston 2020-05-01 275 207.857143 NA NA
Q - Houston 2020-05-02 251 213.142857 NA NA
Q - Houston 2020-05-03 243 223.428571 NA NA
Q - Houston 2020-05-04 182 226.428571 NA NA
Q - Houston 2020-05-05 136 226.000000 NA NA
Q - Houston 2020-05-06 257 231.428571 NA NA
Q - Houston 2020-05-07 181 217.857143 NA NA
Q - Houston 2020-05-08 213 209.000000 NA NA
Q - Houston 2020-05-09 263 210.714286 NA NA
Q - Houston 2020-05-10 224 208.000000 NA NA
Q - Houston 2020-05-11 92 195.142857 NA NA
Q - Houston 2020-05-12 344 224.857143 NA NA
Q - Houston 2020-05-13 310 232.428571 NA NA
Q - Houston 2020-05-14 240 240.857143 NA NA
Q - Houston 2020-05-15 243 245.142857 NA NA
Q - Houston 2020-05-16 332 255.000000 NA NA
Q - Houston 2020-05-17 141 243.142857 NA NA
Q - Houston 2020-05-18 343 279.000000 NA NA
Q - Houston 2020-05-19 214 260.428571 NA NA
Q - Houston 2020-05-20 391 272.000000 NA NA
Q - Houston 2020-05-21 167 261.571429 NA NA
Q - Houston 2020-05-22 243 261.571429 NA NA
Q - Houston 2020-05-23 288 255.285714 NA NA
Q - Houston 2020-05-24 275 274.428571 NA NA
Q - Houston 2020-05-25 171 249.857143 NA NA
Q - Houston 2020-05-26 126 237.285714 NA NA
Q - Houston 2020-05-27 561 261.571429 NA NA
Q - Houston 2020-05-28 463 303.857143 NA NA
Q - Houston 2020-05-29 262 306.571429 NA NA
Q - Houston 2020-05-30 379 319.571429 NA NA
Q - Houston 2020-05-31 752 387.714286 NA NA
Q - Houston 2020-06-01 83 375.142857 NA NA
Q - Houston 2020-06-02 526 432.285714 NA NA
Q - Houston 2020-06-03 603 438.285714 NA NA
Q - Houston 2020-06-04 375 425.714286 NA NA
Q - Houston 2020-06-05 484 457.428571 NA NA
Q - Houston 2020-06-06 446 467.000000 NA NA
Q - Houston 2020-06-07 489 429.428571 NA NA
Q - Houston 2020-06-08 163 440.857143 NA NA
Q - Houston 2020-06-09 416 425.142857 NA NA
Q - Houston 2020-06-10 424 399.571429 NA NA
Q - Houston 2020-06-11 455 411.000000 NA NA
Q - Houston 2020-06-12 385 396.857143 NA NA
Q - Houston 2020-06-13 427 394.142857 NA NA
Q - Houston 2020-06-14 415 383.571429 NA NA
Q - Houston 2020-06-15 221 391.857143 NA NA
Q - Houston 2020-06-16 572 414.142857 NA NA
Q - Houston 2020-06-17 606 440.142857 NA NA
Q - Houston 2020-06-18 708 476.285714 NA NA
Q - Houston 2020-06-19 544 499.000000 NA NA
Q - Houston 2020-06-20 1443 644.142857 NA NA
Q - Houston 2020-06-21 1224 759.714286 NA NA
Q - Houston 2020-06-22 305 771.714286 NA NA
Q - Houston 2020-06-23 2187 1002.428571 NA NA
Q - Houston 2020-06-24 1581 1141.714286 NA NA
Q - Houston 2020-06-25 1575 1265.571429 NA NA
Q - Houston 2020-06-26 1438 1393.285714 NA NA
Q - Houston 2020-06-27 1604 1416.285714 NA NA
Q - Houston 2020-06-28 1019 1387.000000 NA NA
Q - Houston 2020-06-29 142 1363.714286 NA NA
Q - Houston 2020-06-30 1569 1275.428571 NA NA
Q - Houston 2020-07-01 953 1185.714286 NA NA
Q - Houston 2020-07-02 1713 1205.428571 NA NA
Q - Houston 2020-07-03 1554 1222.000000 NA NA
Q - Houston 2020-07-04 1356 1186.571429 NA NA
Q - Houston 2020-07-05 594 1125.857143 NA NA
Q - Houston 2020-07-06 734 1210.428571 NA NA
Q - Houston 2020-07-07 1580 1212.000000 NA NA
Q - Houston 2020-07-08 1785 1330.857143 NA NA
Q - Houston 2020-07-09 996 1228.428571 NA NA
Q - Houston 2020-07-10 1250 1185.000000 NA NA
Q - Houston 2020-07-11 1383 1188.857143 NA NA
Q - Houston 2020-07-12 2137 1409.285714 NA NA
Q - Houston 2020-07-13 1524 1522.142857 NA NA
Q - Houston 2020-07-14 2595 1667.142857 NA NA
Q - Houston 2020-07-15 2253 1734.000000 NA NA
Q - Houston 2020-07-16 2314 1922.285714 NA NA
Q - Houston 2020-07-17 1935 2020.142857 NA NA
Q - Houston 2020-07-18 2205 2137.571429 NA NA
Q - Houston 2020-07-19 1584 2058.571429 NA NA
Q - Houston 2020-07-20 NA 2091.922026 2042.391 2141.453
Q - Houston 2020-07-21 NA 2065.852865 1980.605 2151.101
Q - Houston 2020-07-22 NA 2121.931696 1999.643 2244.221
Q - Houston 2020-07-23 NA 2112.229200 1959.276 2265.182
Q - Houston 2020-07-24 NA 2172.426339 1988.862 2355.990
Q - Houston 2020-07-25 NA 2169.065708 1959.212 2378.920
Q - Houston 2020-07-26 NA 2228.522369 1992.159 2464.886
Q - Houston 2020-07-27 NA 2228.647882 1968.733 2488.562
Q - Houston 2020-07-28 NA 2286.200480 2002.276 2570.125
Q - Houston 2020-07-29 NA 2288.901534 1983.099 2594.704 
#nlme::gapply(TSA.arima.df, FUN = forecast.plot, groups=TSA.daily$Level_Name)
forecast.plot(TSA_Q.arima.df, "TSA Q - Greater Houston Daily Cases", "Daily Cases")
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 138 row(s) containing missing values (geom_path).

## Warning: Removed 138 row(s) containing missing values (geom_path).

######### View Output Table & Graph for Texas ##########
texas.arima.df.1%>%kable(caption="Texas Arima - Daily Cases")%>%kable_styling(full_width=FALSE)
Texas Arima - Daily Cases
Date Cases_Raw Cases_Daily CI_Lower CI_Upper
2020-03-04 0 0.000000 NA NA
2020-03-05 0 0.000000 NA NA
2020-03-06 5 1.666667 NA NA
2020-03-07 0 1.250000 NA NA
2020-03-08 0 1.000000 NA NA
2020-03-09 7 2.000000 NA NA
2020-03-10 3 2.142857 NA NA
2020-03-11 3 2.571429 NA NA
2020-03-12 4 3.142857 NA NA
2020-03-13 0 2.428571 NA NA
2020-03-14 0 2.428571 NA NA
2020-03-15 34 7.285714 NA NA
2020-03-16 0 6.285714 NA NA
2020-03-17 7 6.857143 NA NA
2020-03-18 19 9.142857 NA NA
2020-03-19 26 12.285714 NA NA
2020-03-20 67 21.857143 NA NA
2020-03-21 60 30.428571 NA NA
2020-03-22 28 29.571429 NA NA
2020-03-23 24 33.000000 NA NA
2020-03-24 425 92.714286 NA NA
2020-03-25 263 127.571429 NA NA
2020-03-26 419 183.714286 NA NA
2020-03-27 337 222.285714 NA NA
2020-03-28 317 259.000000 NA NA
2020-03-29 504 327.000000 NA NA
2020-03-30 322 369.571429 NA NA
2020-03-31 392 364.857143 NA NA
2020-04-01 730 431.571429 NA NA
2020-04-02 669 467.285714 NA NA
2020-04-03 659 513.285714 NA NA
2020-04-04 788 580.571429 NA NA
2020-04-05 681 605.857143 NA NA
2020-04-06 484 629.000000 NA NA
2020-04-07 988 714.142857 NA NA
2020-04-08 1092 765.857143 NA NA
2020-04-09 877 795.571429 NA NA
2020-04-10 1441 907.285714 NA NA
2020-04-11 890 921.857143 NA NA
2020-04-12 923 956.428571 NA NA
2020-04-13 422 947.571429 NA NA
2020-04-14 718 909.000000 NA NA
2020-04-15 868 877.000000 NA NA
2020-04-16 963 889.285714 NA NA
2020-04-17 916 814.285714 NA NA
2020-04-18 889 814.142857 NA NA
2020-04-19 663 777.000000 NA NA
2020-04-20 535 793.142857 NA NA
2020-04-21 738 796.000000 NA NA
2020-04-22 873 796.714286 NA NA
2020-04-23 876 784.285714 NA NA
2020-04-24 862 776.571429 NA NA
2020-04-25 968 787.857143 NA NA
2020-04-26 859 815.857143 NA NA
2020-04-27 666 834.571429 NA NA
2020-04-28 874 854.000000 NA NA
2020-04-29 885 855.714286 NA NA
2020-04-30 1033 878.142857 NA NA
2020-05-01 1142 918.142857 NA NA
2020-05-02 1293 964.571429 NA NA
2020-05-03 1026 988.428571 NA NA
2020-05-04 784 1005.285714 NA NA
2020-05-05 1037 1028.571429 NA NA
2020-05-06 1053 1052.571429 NA NA
2020-05-07 1135 1067.142857 NA NA
2020-05-08 1219 1078.142857 NA NA
2020-05-09 1251 1072.142857 NA NA
2020-05-10 1009 1069.714286 NA NA
2020-05-11 1000 1100.571429 NA NA
2020-05-12 1179 1120.857143 NA NA
2020-05-13 1355 1164.000000 NA NA
2020-05-14 1448 1208.714286 NA NA
2020-05-15 1347 1227.000000 NA NA
2020-05-16 1801 1305.571429 NA NA
2020-05-17 785 1273.571429 NA NA
2020-05-18 909 1260.571429 NA NA
2020-05-19 1219 1266.285714 NA NA
2020-05-20 1411 1274.285714 NA NA
2020-05-21 945 1202.428571 NA NA
2020-05-22 1241 1187.285714 NA NA
2020-05-23 1060 1081.428571 NA NA
2020-05-24 839 1089.142857 NA NA
2020-05-25 623 1048.285714 NA NA
2020-05-26 620 962.714286 NA NA
2020-05-27 1361 955.571429 NA NA
2020-05-28 1855 1085.571429 NA NA
2020-05-29 1230 1084.000000 NA NA
2020-05-30 1333 1123.000000 NA NA
2020-05-31 1949 1281.571429 NA NA
2020-06-01 593 1277.285714 NA NA
2020-06-02 1688 1429.857143 NA NA
2020-06-03 1703 1478.714286 NA NA
2020-06-04 1650 1449.428571 NA NA
2020-06-05 1693 1515.571429 NA NA
2020-06-06 1942 1602.571429 NA NA
2020-06-07 1426 1527.857143 NA NA
2020-06-08 638 1534.285714 NA NA
2020-06-09 1853 1557.857143 NA NA
2020-06-10 2509 1673.000000 NA NA
2020-06-11 1826 1698.142857 NA NA
2020-06-12 2097 1755.857143 NA NA
2020-06-13 2331 1811.428571 NA NA
2020-06-14 1843 1871.000000 NA NA
2020-06-15 1254 1959.000000 NA NA
2020-06-16 4098 2279.714286 NA NA
2020-06-17 3140 2369.857143 NA NA
2020-06-18 3516 2611.285714 NA NA
2020-06-19 3454 2805.142857 NA NA
2020-06-20 4430 3105.000000 NA NA
2020-06-21 3866 3394.000000 NA NA
2020-06-22 3280 3683.428571 NA NA
2020-06-23 5522 3886.857143 NA NA
2020-06-24 5551 4231.285714 NA NA
2020-06-25 5996 4585.571429 NA NA
2020-06-26 5707 4907.428571 NA NA
2020-06-27 5742 5094.857143 NA NA
2020-06-28 5357 5307.857143 NA NA
2020-06-29 4288 5451.857143 NA NA
2020-06-30 6975 5659.428571 NA NA
2020-07-01 8076 6020.142857 NA NA
2020-07-02 7915 6294.285714 NA NA
2020-07-03 7555 6558.285714 NA NA
2020-07-04 8258 6917.714286 NA NA
2020-07-05 3449 6645.142857 NA NA
2020-07-06 5319 6792.428571 NA NA
2020-07-07 10028 7228.571429 NA NA
2020-07-08 9979 7500.428571 NA NA
2020-07-09 9782 7767.142857 NA NA
2020-07-10 9765 8082.857143 NA NA
2020-07-11 10351 8381.857143 NA NA
2020-07-12 8196 9060.000000 NA NA
2020-07-13 5655 9108.000000 NA NA
2020-07-14 10745 9210.428571 NA NA
2020-07-15 9550 9149.142857 NA NA
2020-07-16 10291 9221.857143 NA NA
2020-07-17 14916 9957.714286 NA NA
2020-07-18 10344 9956.714286 NA NA
2020-07-19 7322 9831.857143 NA NA
2020-07-20 NA 10001.122962 9884.249 10118.00
2020-07-21 NA 10170.388780 9991.307 10349.47
2020-07-22 NA 10339.654599 10103.139 10576.17
2020-07-23 NA 10508.920418 10215.704 10802.14
2020-07-24 NA 10678.186237 10327.634 11028.74
2020-07-25 NA 10847.452055 10438.332 11256.57
2020-07-26 NA 11016.717874 10547.504 11485.93
2020-07-27 NA 11185.983693 10655.001 11716.97
2020-07-28 NA 11355.249512 10760.751 11949.75
2020-07-29 NA 11524.515330 10864.722 12184.31 
forecast.plot(texas.arima.df.1, "Texas Daily Cases", "Daily Cases")
## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 138 row(s) containing missing values (geom_path).

## Warning: Removed 138 row(s) containing missing values (geom_path).

export

write.csv(texas.arima.df.1, 'statistical-output/time-series/state_case_timeseries.csv', row.names = FALSE)
write.csv(TSA.arima.df.1, 'statistical-output/time-series/tsa_case_timeseries.csv', row.names = FALSE)
write.csv(county.arima.df.1, 'statistical-output/time-series/county_case_timeseries.csv', row.names = FALSE)
write.csv(metro.arima.df.1, 'statistical-output/time-series/metro_case_timeseries.csv', row.names = FALSE)
write.csv(PHR.arima.df, 'statistical-output/time-series/phr_case_timeseries.csv', row.names = FALSE)

# write.csv(texas.arima.df.2, 'statistical-output/time-series/state_timeseries_06_03.csv', row.names = FALSE)
# write.csv(TSA.arima.df.2, 'statistical-output/time-series/tsa_timeseries_06_03.csv', row.names = FALSE)
# write.csv(county.arima.df.2, 'statistical-output/time-series/county_timeseries_06_03.csv', row.names = FALSE)
# write.csv(metro.arima.df.2, 'statistical-output/time-series/metro_timeseries_06_03.csv', row.names = FALSE)

grouping

colnames(county.arima.df.1)[1] <- 'Level'
colnames(metro.arima.df.1)[1] <- 'Level'
colnames(TSA.arima.df.1)[1] <- 'Level'
texas.arima.df.1$Level <- 'Texas'
colnames(PHR.arima.df)[1] <- 'Level'

county.arima.df.1$Level_Type = 'County'
metro.arima.df.1$Level_Type = 'Metro'
TSA.arima.df.1$Level_Type = 'TSA'
PHR.arima.df$Level_Type = 'PHR'
texas.arima.df.1$Level_Type = 'State'


combined.arima.df.1 <- rbind(county.arima.df.1, metro.arima.df.1, TSA.arima.df.1, texas.arima.df.1, PHR.arima.df)
write.csv(combined.arima.df.1, 'statistical-output/time-series/stacked_case_timeseries.csv',
          row.names = FALSE)

write.csv(combined.arima.df.1, 'tableau/stacked_case_timeseries.csv',
          row.names = FALSE)
# 
# combined.arima.df.2 <- rbind(county.arima.df.2, metro.arima.df.2, TSA.arima.df.2, texas.arima.df.2)
# write.csv(combined.arima.df.2, 'statistical-output/time-series/stacked_timeseries_06_03.csv',
#           row.names = FALSE)
# 
# write.csv(combined.arima.df.2, 'tableau/stacked_timeseries_06_03.csv',
#           row.names = FALSE)

Hospitalization Time Series

# Compute forecast (UPDATE PREDICTION PERIOD [days] AS NEEDED)
covid.arima.forecast<-function(mydata, prediction.period = 10, mindate)
{
  maxdate <- max(mydata$Date)
  # mindate <- as.Date('2020-05-01')
  pred_start_label = format(mindate, format = '%m_%d')
  
  mydata = subset(mydata, Date >= mindate)
  model.length <- as.numeric(length(mydata$Date) + prediction.period)

  if(max(mydata$Hospitalizations_Total>=100, na.rm = TRUE))
  {

    
    # arima requires cases to be a timeseries vector
    #convert daily cases to time series
    my.timeseries<-ts(mydata$Hospitalizations_Total)
    #my.timeseries <- rollmeanr(my.timeseries, k=7, na.pad=TRUE, align = 'right')
    
    #load package(pracma)
    library(pracma)
    my.timeseries<-movavg(my.timeseries,7,"s")
    #d=0 restricts first differencing to 0 so that daily cases aren't differenced
    
    arima.fit <- forecast::auto.arima(my.timeseries)
  
    # obtain diagnostic plots for ideal arima (p,d,q) selection 
    acf <- ggAcf(my.timeseries, lag.max=30)
    pacf <- ggPacf(my.timeseries, lag.max=30)
    ts.diagnostics <- grid.arrange(acf, pacf, nrow=2)
    ggsave(paste0('statistical-output/time-series/diagnostics/',
                  mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'),
           plot = ts.diagnostics)

    
    
    # save parameters from arima autofit
    p<- arima.fit$arma[1]          # autoregressive order 
    q<- arima.fit$arma[2]          # moving average order 
    d<-arima.fit$arma[6]           # differencing order from model
  
    # 10 day forecast, CI for lower and upper has confidence level 70% set by level =c(70,70)
    arima.forecast <- forecast::forecast(arima.fit, h = prediction.period, level=c(70,70))

    #return a dataframe of the arima model(Daily cases by date)
    arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
                            Cases_Raw = c(mydata$Hospitalizations_Total, rep(NA, times = prediction.period)),
                            Hospitalizations_Total = c(my.timeseries, arima.forecast[['mean']]),
                            CI_Lower = c(rep(NA, times = length(my.timeseries)),
                                         arima.forecast[['lower']][, 2]),
                            CI_Upper = c(rep(NA, times = length(my.timeseries)),
                                         arima.forecast[['upper']][, 2]))
                            # Order_AutoReg = c(rep(p, times = model.length)),
                            # Order_Moving_Avg = c(rep(q, times = model.length)),
                            # Differencing = c(rep(d, times = model.length)))
    
      # save prediction plot for preliminary review
      arima.plot <- ggplot(arima.out, aes(x=Date, y = Hospitalizations_Total))+
                    geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), fill = "red", alpha = 0.5, size = 0.1)+
                    geom_line(color = "black", size = 1)+
                    labs(y = 'Daily Cases (7-Day Moving Average)', x = 'Date',
                         title = mydata$Level_Name[1]) +
                    scale_x_date(limits = c(mindate, maxdate + prediction.period),
                                 date_labels = '%b-%d', date_breaks = '1 week') + 
                    ggpubr::theme_pubr() +
                    theme(axis.text.x = element_text(angle = -90))
      
      ggsave(plot = arima.plot, paste0('statistical-output/time-series/plots/',
                                       mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'))    
    
    } else {
    # insufficient data catch: return NA values for predictions 
    arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
                            Cases_Raw = c(mydata$Hospitalizations_Total, rep(NA, times = prediction.period)),
                            Hospitalizations_Total = rep(NA, times = model.length),
                            CI_Lower = rep(NA, times = model.length),
                            CI_Upper =  rep(NA, times = model.length))
                            # Order_AutoReg = c(rep(NA, times = model.length)),
                            # Order_Moving_Avg = c(rep(NA, times = model.length)),
                            # Differencing = c(rep(NA, times= model.length)))
    }
  
  #replace CI lower limit with 0 if negative
  arima.out$CI_Lower <- ifelse(arima.out$CI_Lower>=0,arima.out$CI_Lower, 0)
  
  
  return(arima.out)
}

TSA - Hospitalization

# All data
TSA.arima.output.1 <- nlme::gapply(TSA.Daily,
                                   FUN = covid.arima.forecast,
                                   groups = TSA.Daily$Level_Name,
                                   mindate = as.Date('2020-03-04'))
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 row(s) containing missing values (geom_path).

TSA.arima.df.1 <- data.table::rbindlist(TSA.arima.output.1, idcol=TRUE)

colnames(TSA.arima.df.1)[1] <- 'TSA'
TSA.arima.df.1 <- merge(TSA.arima.df.1, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = c('TSA'))

# combine TSA 
TSA.arima.df.1$TSA <- paste0(TSA.arima.df.1$TSA, ' - ', TSA.arima.df.1$TSA_Name)
TSA.arima.df.1$TSA_Name <- NULL

State - hospitalization

texas.arima.df.1 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-03-04'))
## Saving 7 x 5 in image
## Saving 7 x 5 in image

# texas.arima.df.2 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-06-03'))

Hospitalization Write to CSV

tsa.hosp.arima<-write.csv(TSA.arima.df.1, "statistical-output/time-series/tsa_hosp_timeseries.csv", row.names=FALSE)
texas.hosp.arima<-write.csv(texas.arima.df.1, "statistical-output/time-series/state_hosp_timeseries.csv", row.names=FALSE)

Stacking

TSA_hosp_out = TSA.arima.df.1
colnames(TSA_hosp_out)[1] = 'Level'
texas.arima.df.1$Level = 'Texas'


TSA_hosp_out$Level_Type = 'TSA'
texas.arima.df.1$Level_Type = 'State'

combined.arima.df.1 <- rbind(TSA_hosp_out, texas.arima.df.1)
write.csv(combined.arima.df.1, 'statistical-output/time-series/stacked_hosp_timeseries.csv',
          row.names = FALSE)

write.csv(combined.arima.df.1, 'tableau/stacked_hosp_timeseries.csv',
          row.names = FALSE)

Hospitalization Forecast Plots

#Create plot function, forecast.data is the dataframe containing true data
# and the forecast data. days.forecast is the number of days that are forecast
# in the dataframe (here we are forecasting 10 days, so it should be 10 unless we
# change this at some point down the road)
forecast.plot<-function(forecast.data, plot.title, y.label){
    mindate<-as.POSIXct("2020-03-04")
    maxdate<-as.POSIXct(max(forecast.data$Date))
  ts.plot<-ggplot(aes(x=as.POSIXct(Date), y=Hospitalizations_Total), data=forecast.data)+
    geom_ribbon(aes(ymin=CI_Lower, ymax=CI_Upper), fill="grey50", size=0.1)+
    geom_line(color="blue", size=1)+
    geom_line(aes(x=as.POSIXct(Date), y=CI_Lower), color="grey", size=0.1)+
    geom_line(aes(x=as.POSIXct(Date),y=CI_Upper), color="grey", size=0.1)+
    scale_x_datetime(limits = c(mindate, maxdate))+
    xlab("Date")+
    ylab(y.label)+
    #Can use the following title if we are running using nlme for data frame
    #ggtitle(paste(forecast.data$.id,": TS Daily Cases + 10 Day Forcast",sep=""))
    ggtitle(plot.title)
ts.plot
}
######### View Output Table & Graph for TSA Q ##########
library(kableExtra)
#subset TSA Q arima data
TSA_Q.arima.df<-TSA.arima.df.1%>%subset(TSA.arima.df.1$TSA == "Q - Houston")
#apply the plot function to TSA Q arima data frame
TSA_Q.arima.df%>%kable(caption="Greater Houston ARIMA - Daily Hosp")%>%kable_styling(full_width=FALSE)
Greater Houston ARIMA - Daily Hosp
TSA Date Cases_Raw Hospitalizations_Total CI_Lower CI_Upper
Q - Houston 2020-03-04 NA NA NA NA
Q - Houston 2020-03-05 NA NA NA NA
Q - Houston 2020-03-06 NA NA NA NA
Q - Houston 2020-03-07 NA NA NA NA
Q - Houston 2020-03-08 NA NA NA NA
Q - Houston 2020-03-09 NA NA NA NA
Q - Houston 2020-03-10 NA NA NA NA
Q - Houston 2020-03-11 NA NA NA NA
Q - Houston 2020-03-12 NA NA NA NA
Q - Houston 2020-03-13 NA NA NA NA
Q - Houston 2020-03-14 NA NA NA NA
Q - Houston 2020-03-15 NA NA NA NA
Q - Houston 2020-03-16 NA NA NA NA
Q - Houston 2020-03-17 NA NA NA NA
Q - Houston 2020-03-18 NA NA NA NA
Q - Houston 2020-03-19 NA NA NA NA
Q - Houston 2020-03-20 NA NA NA NA
Q - Houston 2020-03-21 NA NA NA NA
Q - Houston 2020-03-22 NA NA NA NA
Q - Houston 2020-03-23 NA NA NA NA
Q - Houston 2020-03-24 NA NA NA NA
Q - Houston 2020-03-25 NA NA NA NA
Q - Houston 2020-03-26 NA NA NA NA
Q - Houston 2020-03-27 NA NA NA NA
Q - Houston 2020-03-28 NA NA NA NA
Q - Houston 2020-03-29 NA NA NA NA
Q - Houston 2020-03-30 NA NA NA NA
Q - Houston 2020-03-31 NA NA NA NA
Q - Houston 2020-04-01 NA NA NA NA
Q - Houston 2020-04-02 NA NA NA NA
Q - Houston 2020-04-03 NA NA NA NA
Q - Houston 2020-04-04 NA NA NA NA
Q - Houston 2020-04-05 NA NA NA NA
Q - Houston 2020-04-06 NA NA NA NA
Q - Houston 2020-04-07 NA NA NA NA
Q - Houston 2020-04-08 NA NA NA NA
Q - Houston 2020-04-09 NA NA NA NA
Q - Houston 2020-04-10 NA NA NA NA
Q - Houston 2020-04-11 NA NA NA NA
Q - Houston 2020-04-12 516 NA NA NA
Q - Houston 2020-04-13 319 NA NA NA
Q - Houston 2020-04-14 527 NA NA NA
Q - Houston 2020-04-15 584 NA NA NA
Q - Houston 2020-04-16 511 NA NA NA
Q - Houston 2020-04-17 500 NA NA NA
Q - Houston 2020-04-18 363 474.2857 NA NA
Q - Houston 2020-04-19 535 477.0000 NA NA
Q - Houston 2020-04-20 516 505.1429 NA NA
Q - Houston 2020-04-21 516 503.5714 NA NA
Q - Houston 2020-04-22 505 492.2857 NA NA
Q - Houston 2020-04-23 501 490.8571 NA NA
Q - Houston 2020-04-24 483 488.4286 NA NA
Q - Houston 2020-04-25 465 503.0000 NA NA
Q - Houston 2020-04-26 464 492.8571 NA NA
Q - Houston 2020-04-27 468 486.0000 NA NA
Q - Houston 2020-04-28 485 481.5714 NA NA
Q - Houston 2020-04-29 436 471.7143 NA NA
Q - Houston 2020-04-30 442 463.2857 NA NA
Q - Houston 2020-05-01 475 462.1429 NA NA
Q - Houston 2020-05-02 434 457.7143 NA NA
Q - Houston 2020-05-03 329 438.4286 NA NA
Q - Houston 2020-05-04 418 431.2857 NA NA
Q - Houston 2020-05-05 544 439.7143 NA NA
Q - Houston 2020-05-06 436 439.7143 NA NA
Q - Houston 2020-05-07 453 441.2857 NA NA
Q - Houston 2020-05-08 430 434.8571 NA NA
Q - Houston 2020-05-09 446 436.5714 NA NA
Q - Houston 2020-05-10 446 453.2857 NA NA
Q - Houston 2020-05-11 402 451.0000 NA NA
Q - Houston 2020-05-12 436 435.5714 NA NA
Q - Houston 2020-05-13 435 435.4286 NA NA
Q - Houston 2020-05-14 411 429.4286 NA NA
Q - Houston 2020-05-15 401 425.2857 NA NA
Q - Houston 2020-05-16 414 420.7143 NA NA
Q - Houston 2020-05-17 380 411.2857 NA NA
Q - Houston 2020-05-18 397 410.5714 NA NA
Q - Houston 2020-05-19 433 410.1429 NA NA
Q - Houston 2020-05-20 443 411.2857 NA NA
Q - Houston 2020-05-21 444 416.0000 NA NA
Q - Houston 2020-05-22 406 416.7143 NA NA
Q - Houston 2020-05-23 387 412.8571 NA NA
Q - Houston 2020-05-24 407 416.7143 NA NA
Q - Houston 2020-05-25 392 416.0000 NA NA
Q - Houston 2020-05-26 401 411.4286 NA NA
Q - Houston 2020-05-27 414 407.2857 NA NA
Q - Houston 2020-05-28 435 406.0000 NA NA
Q - Houston 2020-05-29 407 406.1429 NA NA
Q - Houston 2020-05-30 460 416.5714 NA NA
Q - Houston 2020-05-31 454 423.2857 NA NA
Q - Houston 2020-06-01 482 436.1429 NA NA
Q - Houston 2020-06-02 506 451.1429 NA NA
Q - Houston 2020-06-03 501 463.5714 NA NA
Q - Houston 2020-06-04 526 476.5714 NA NA
Q - Houston 2020-06-05 555 497.7143 NA NA
Q - Houston 2020-06-06 522 506.5714 NA NA
Q - Houston 2020-06-07 565 522.4286 NA NA
Q - Houston 2020-06-08 614 541.2857 NA NA
Q - Houston 2020-06-09 622 557.8571 NA NA
Q - Houston 2020-06-10 690 584.8571 NA NA
Q - Houston 2020-06-11 642 601.4286 NA NA
Q - Houston 2020-06-12 663 616.8571 NA NA
Q - Houston 2020-06-13 681 639.5714 NA NA
Q - Houston 2020-06-14 696 658.2857 NA NA
Q - Houston 2020-06-15 730 674.8571 NA NA
Q - Houston 2020-06-16 795 699.5714 NA NA
Q - Houston 2020-06-17 839 720.8571 NA NA
Q - Houston 2020-06-18 940 763.4286 NA NA
Q - Houston 2020-06-19 1031 816.0000 NA NA
Q - Houston 2020-06-20 1040 867.2857 NA NA
Q - Houston 2020-06-21 1062 919.5714 NA NA
Q - Houston 2020-06-22 1185 984.5714 NA NA
Q - Houston 2020-06-23 1293 1055.7143 NA NA
Q - Houston 2020-06-24 1343 1127.7143 NA NA
Q - Houston 2020-06-25 1496 1207.1429 NA NA
Q - Houston 2020-06-26 1606 1289.2857 NA NA
Q - Houston 2020-06-27 1691 1382.2857 NA NA
Q - Houston 2020-06-28 1642 1465.1429 NA NA
Q - Houston 2020-06-29 1767 1548.2857 NA NA
Q - Houston 2020-06-30 2059 1657.7143 NA NA
Q - Houston 2020-07-01 2139 1771.4286 NA NA
Q - Houston 2020-07-02 2250 1879.1429 NA NA
Q - Houston 2020-07-03 2257 1972.1429 NA NA
Q - Houston 2020-07-04 2369 2069.0000 NA NA
Q - Houston 2020-07-05 2442 2183.2857 NA NA
Q - Houston 2020-07-06 2656 2310.2857 NA NA
Q - Houston 2020-07-07 2699 2401.7143 NA NA
Q - Houston 2020-07-08 2700 2481.8571 NA NA
Q - Houston 2020-07-09 2742 2552.1429 NA NA
Q - Houston 2020-07-10 2826 2633.4286 NA NA
Q - Houston 2020-07-11 2842 2701.0000 NA NA
Q - Houston 2020-07-12 2863 2761.1429 NA NA
Q - Houston 2020-07-13 2842 2787.7143 NA NA
Q - Houston 2020-07-14 2914 2818.4286 NA NA
Q - Houston 2020-07-15 2884 2844.7143 NA NA
Q - Houston 2020-07-16 2815 2855.1429 NA NA
Q - Houston 2020-07-17 2763 2846.1429 NA NA
Q - Houston 2020-07-18 2785 2838.0000 NA NA
Q - Houston 2020-07-19 2742 2820.7143 NA NA
Q - Houston 2020-07-20 NA 2803.4286 2791.592 2815.265
Q - Houston 2020-07-21 NA 2786.1429 2759.675 2812.610
Q - Houston 2020-07-22 NA 2768.8571 2724.569 2813.146
Q - Houston 2020-07-23 NA 2751.5714 2686.740 2816.403
Q - Houston 2020-07-24 NA 2734.2857 2646.503 2822.068
Q - Houston 2020-07-25 NA 2717.0000 2604.086 2829.914
Q - Houston 2020-07-26 NA 2699.7143 2559.661 2839.767
Q - Houston 2020-07-27 NA 2682.4286 2513.368 2851.489
Q - Houston 2020-07-28 NA 2665.1429 2465.318 2864.968
Q - Houston 2020-07-29 NA 2647.8571 2415.606 2880.109 
#nlme::gapply(TSA.arima.df, FUN = forecast.plot, groups=TSA.daily$Level_Name)
forecast.plot(TSA_Q.arima.df, "TSA Q - Greater Houston Daily Hosp", "Daily Hosp")
## Warning: Removed 45 row(s) containing missing values (geom_path).
## Warning: Removed 138 row(s) containing missing values (geom_path).

## Warning: Removed 138 row(s) containing missing values (geom_path).

######### View Output Table & Graph for Texas ##########
texas.arima.df.1%>%kable(caption="Texas Arima - Daily Hosp")%>%kable_styling(full_width=FALSE)
Texas Arima - Daily Hosp
Date Cases_Raw Hospitalizations_Total CI_Lower CI_Upper Level Level_Type
2020-03-04 0 0.0000 NA NA Texas State
2020-03-05 0 0.0000 NA NA Texas State
2020-03-06 0 0.0000 NA NA Texas State
2020-03-07 0 0.0000 NA NA Texas State
2020-03-08 0 0.0000 NA NA Texas State
2020-03-09 0 0.0000 NA NA Texas State
2020-03-10 0 0.0000 NA NA Texas State
2020-03-11 0 0.0000 NA NA Texas State
2020-03-12 0 0.0000 NA NA Texas State
2020-03-13 0 0.0000 NA NA Texas State
2020-03-14 0 0.0000 NA NA Texas State
2020-03-15 0 0.0000 NA NA Texas State
2020-03-16 0 0.0000 NA NA Texas State
2020-03-17 0 0.0000 NA NA Texas State
2020-03-18 0 0.0000 NA NA Texas State
2020-03-19 0 0.0000 NA NA Texas State
2020-03-20 0 0.0000 NA NA Texas State
2020-03-21 0 0.0000 NA NA Texas State
2020-03-22 0 0.0000 NA NA Texas State
2020-03-23 0 0.0000 NA NA Texas State
2020-03-24 0 0.0000 NA NA Texas State
2020-03-25 0 0.0000 NA NA Texas State
2020-03-26 0 0.0000 NA NA Texas State
2020-03-27 0 0.0000 NA NA Texas State
2020-03-28 0 0.0000 NA NA Texas State
2020-03-29 0 0.0000 NA NA Texas State
2020-03-30 0 0.0000 NA NA Texas State
2020-03-31 0 0.0000 NA NA Texas State
2020-04-01 0 0.0000 NA NA Texas State
2020-04-02 0 0.0000 NA NA Texas State
2020-04-03 0 0.0000 NA NA Texas State
2020-04-04 0 0.0000 NA NA Texas State
2020-04-05 0 0.0000 NA NA Texas State
2020-04-06 0 0.0000 NA NA Texas State
2020-04-07 0 0.0000 NA NA Texas State
2020-04-08 0 0.0000 NA NA Texas State
2020-04-09 0 0.0000 NA NA Texas State
2020-04-10 0 0.0000 NA NA Texas State
2020-04-11 0 0.0000 NA NA Texas State
2020-04-12 1338 191.1429 NA NA Texas State
2020-04-13 1176 359.1429 NA NA Texas State
2020-04-14 1409 560.4286 NA NA Texas State
2020-04-15 1568 784.4286 NA NA Texas State
2020-04-16 1459 992.8571 NA NA Texas State
2020-04-17 1522 1210.2857 NA NA Texas State
2020-04-18 1321 1399.0000 NA NA Texas State
2020-04-19 1471 1418.0000 NA NA Texas State
2020-04-20 1414 1452.0000 NA NA Texas State
2020-04-21 1497 1464.5714 NA NA Texas State
2020-04-22 1678 1480.2857 NA NA Texas State
2020-04-23 1649 1507.4286 NA NA Texas State
2020-04-24 1664 1527.7143 NA NA Texas State
2020-04-25 1587 1565.7143 NA NA Texas State
2020-04-26 1542 1575.8571 NA NA Texas State
2020-04-27 1563 1597.1429 NA NA Texas State
2020-04-28 1682 1623.5714 NA NA Texas State
2020-04-29 1702 1627.0000 NA NA Texas State
2020-04-30 1686 1632.2857 NA NA Texas State
2020-05-01 1778 1648.5714 NA NA Texas State
2020-05-02 1725 1668.2857 NA NA Texas State
2020-05-03 1540 1668.0000 NA NA Texas State
2020-05-04 1533 1663.7143 NA NA Texas State
2020-05-05 1888 1693.1429 NA NA Texas State
2020-05-06 1812 1708.8571 NA NA Texas State
2020-05-07 1750 1718.0000 NA NA Texas State
2020-05-08 1734 1711.7143 NA NA Texas State
2020-05-09 1735 1713.1429 NA NA Texas State
2020-05-10 1626 1725.4286 NA NA Texas State
2020-05-11 1525 1724.2857 NA NA Texas State
2020-05-12 1725 1701.0000 NA NA Texas State
2020-05-13 1676 1681.5714 NA NA Texas State
2020-05-14 1648 1667.0000 NA NA Texas State
2020-05-15 1716 1664.4286 NA NA Texas State
2020-05-16 1791 1672.4286 NA NA Texas State
2020-05-17 1512 1656.1429 NA NA Texas State
2020-05-18 1551 1659.8571 NA NA Texas State
2020-05-19 1732 1660.8571 NA NA Texas State
2020-05-20 1791 1677.2857 NA NA Texas State
2020-05-21 1680 1681.8571 NA NA Texas State
2020-05-22 1578 1662.1429 NA NA Texas State
2020-05-23 1688 1647.4286 NA NA Texas State
2020-05-24 1572 1656.0000 NA NA Texas State
2020-05-25 1511 1650.2857 NA NA Texas State
2020-05-26 1534 1622.0000 NA NA Texas State
2020-05-27 1645 1601.1429 NA NA Texas State
2020-05-28 1692 1602.8571 NA NA Texas State
2020-05-29 1622 1609.1429 NA NA Texas State
2020-05-30 1752 1618.2857 NA NA Texas State
2020-05-31 1684 1634.2857 NA NA Texas State
2020-06-01 1756 1669.2857 NA NA Texas State
2020-06-02 1773 1703.4286 NA NA Texas State
2020-06-03 1805 1726.2857 NA NA Texas State
2020-06-04 1796 1741.1429 NA NA Texas State
2020-06-05 1855 1774.4286 NA NA Texas State
2020-06-06 1822 1784.4286 NA NA Texas State
2020-06-07 1878 1812.1429 NA NA Texas State
2020-06-08 1935 1837.7143 NA NA Texas State
2020-06-09 2056 1878.1429 NA NA Texas State
2020-06-10 2153 1927.8571 NA NA Texas State
2020-06-11 2059 1965.4286 NA NA Texas State
2020-06-12 2166 2009.8571 NA NA Texas State
2020-06-13 2242 2069.8571 NA NA Texas State
2020-06-14 2287 2128.2857 NA NA Texas State
2020-06-15 2326 2184.1429 NA NA Texas State
2020-06-16 2518 2250.1429 NA NA Texas State
2020-06-17 2793 2341.5714 NA NA Texas State
2020-06-18 2947 2468.4286 NA NA Texas State
2020-06-19 3148 2608.7143 NA NA Texas State
2020-06-20 3247 2752.2857 NA NA Texas State
2020-06-21 3409 2912.5714 NA NA Texas State
2020-06-22 3711 3110.4286 NA NA Texas State
2020-06-23 4092 3335.2857 NA NA Texas State
2020-06-24 4389 3563.2857 NA NA Texas State
2020-06-25 4739 3819.2857 NA NA Texas State
2020-06-26 5102 4098.4286 NA NA Texas State
2020-06-27 5523 4423.5714 NA NA Texas State
2020-06-28 5497 4721.8571 NA NA Texas State
2020-06-29 5913 5036.4286 NA NA Texas State
2020-06-30 6533 5385.1429 NA NA Texas State
2020-07-01 6904 5744.4286 NA NA Texas State
2020-07-02 7382 6122.0000 NA NA Texas State
2020-07-03 7652 6486.2857 NA NA Texas State
2020-07-04 7890 6824.4286 NA NA Texas State
2020-07-05 8181 7207.8571 NA NA Texas State
2020-07-06 8698 7605.7143 NA NA Texas State
2020-07-07 9286 7999.0000 NA NA Texas State
2020-07-08 9610 8385.5714 NA NA Texas State
2020-07-09 9689 8715.1429 NA NA Texas State
2020-07-10 10002 9050.8571 NA NA Texas State
2020-07-11 10083 9364.1429 NA NA Texas State
2020-07-12 10410 9682.5714 NA NA Texas State
2020-07-13 10405 9926.4286 NA NA Texas State
2020-07-14 10569 10109.7143 NA NA Texas State
2020-07-15 10471 10232.7143 NA NA Texas State
2020-07-16 10457 10342.4286 NA NA Texas State
2020-07-17 10632 10432.4286 NA NA Texas State
2020-07-18 10658 10514.5714 NA NA Texas State
2020-07-19 10592 10540.5714 NA NA Texas State
2020-07-20 NA 10566.5714 10536.40 10596.75 Texas State
2020-07-21 NA 10592.5714 10525.10 10660.05 Texas State
2020-07-22 NA 10618.5714 10505.66 10731.48 Texas State
2020-07-23 NA 10644.5714 10479.29 10809.85 Texas State
2020-07-24 NA 10670.5714 10446.78 10894.36 Texas State
2020-07-25 NA 10696.5714 10408.71 10984.43 Texas State
2020-07-26 NA 10722.5714 10365.53 11079.62 Texas State
2020-07-27 NA 10748.5714 10317.58 11179.57 Texas State
2020-07-28 NA 10774.5714 10265.15 11284.00 Texas State
2020-07-29 NA 10800.5714 10208.48 11392.66 Texas State 
forecast.plot(texas.arima.df.1, "Texas Daily Hospitalizations", "Daily Hosp")
## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 138 row(s) containing missing values (geom_path).

## Warning: Removed 138 row(s) containing missing values (geom_path).

STANDARD STATISTICAL TESTS

CASE RATIOS

TSA

#TODO: update all varnames
cumcases.TSA <- TSA.Daily %>% dplyr::select(Date,Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.TSA$Date <- ymd(cumcases.TSA$Date)

latestdate <- max(cumcases.TSA$Date)
earliestdate <- latestdate - 14
middate <- latestdate - 7

week2 <- subset(cumcases.TSA, Date==latestdate, select=Cases_Cumulative) - 
  subset(cumcases.TSA, Date==middate, select=Cases_Cumulative)

week1 <- subset(cumcases.TSA, Date==middate, select=Cases_Cumulative) - 
  subset(cumcases.TSA, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
#current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA
tsa_current_ratio_dat <- data.frame(TSA=subset(cumcases.TSA, Date==latestdate, select=Level_Name)[,1],
                                    current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

# add TSA name
tsa_current_ratio_dat <- merge(tsa_current_ratio_dat, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = 'TSA')

# combine TSA name
tsa_current_ratio_dat$TSA <- paste0(tsa_current_ratio_dat$TSA, ' - ', tsa_current_ratio_dat$TSA_Name)
tsa_current_ratio_dat$TSA_Name <- NULL

PHR

cumcases.PHR <- PHR.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.PHR$Date <- ymd(cumcases.PHR$Date)

latestdate <- max(cumcases.PHR$Date)
earliestdate <- latestdate - 14
middate <- latestdate - 7

week2 <- subset(cumcases.PHR, Date==latestdate, select=Cases_Cumulative) - 
  subset(cumcases.PHR, Date==middate, select=Cases_Cumulative)

week1 <- subset(cumcases.PHR, Date==middate, select=Cases_Cumulative) - 
  subset(cumcases.PHR, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
#current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA
phr_current_ratio_dat <- data.frame(Level_Name=subset(cumcases.PHR, Date==latestdate, select=Level_Name)[,1],
                                    current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

# add PHR name
phr_current_ratio_dat <- merge(phr_current_ratio_dat, unique(PHR.Daily[, c('PHR', 'Level_Name')]), by = 'Level_Name')

# combine PHR name
phr_current_ratio_dat$Level_Name <- paste0(phr_current_ratio_dat$PHR, ' - ', phr_current_ratio_dat$Level_Name)
phr_current_ratio_dat$PHR <- NULL

colnames(phr_current_ratio_dat)[1] <- 'PHR'

County

# declare vals
cumcases.county <- county.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.county$Date <- ymd(cumcases.county$Date)

latestdate <- max(cumcases.county$Date, na.rm=T)
earliestdate <- latestdate - 14
middate <- latestdate-7

# calc cumulative case ratios from 2 weeks ago and last week
week2 <- subset(cumcases.county, Date==latestdate, select=Cases_Cumulative) - 
         subset(cumcases.county, Date==middate, select=Cases_Cumulative)

week1 <- subset(cumcases.county, Date==middate, select=Cases_Cumulative) - 
         subset(cumcases.county, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]

# add contingencies 
current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA

# output
# TODO: investigate select vs subset(...)[, 'County']
county_current_ratio_dat <- 
  data.frame(County=subset(cumcases.county, Date==latestdate, select=Level_Name)[,1],
             current_ratio=current_ratio,
             current_ratio_cat = cut(current_ratio,
                                     breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

Metro

cumcases.metro <- metro.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.metro$Date <- ymd(cumcases.metro$Date)

latestdate <- max(cumcases.metro$Date, na.rm=T)
earliestdate <- latestdate - 14
middate <- latestdate-7

week2 <- subset(cumcases.metro, Date==latestdate, select=Cases_Cumulative) - 
         subset(cumcases.metro, Date==middate, select=Cases_Cumulative)
week1 <- subset(cumcases.metro, Date==middate, select=Cases_Cumulative) - 
         subset(cumcases.metro, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
current_ratio[current_ratio<0] <- NA

metro_current_ratio_dat <-
  data.frame(Metro_Area=subset(cumcases.metro, Date==latestdate, select=Level_Name)[,1],
             current_ratio=current_ratio,
             current_ratio_cat = cut(current_ratio,
                                     breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

State

cumcases.state <- state.Daily %>% dplyr::select(Date, Cases_Cumulative)
cumcases.state$Date <- ymd(cumcases.state$Date)

latestdate <- max(cumcases.state$Date, na.rm = TRUE)
earliestdate <- latestdate - 14
middate <- latestdate-7

week2 <- subset(cumcases.state, Date == latestdate, select = Cases_Cumulative) - 
         subset(cumcases.state, Date == middate, select = Cases_Cumulative)

week1 <- subset(cumcases.state, Date == middate, select = Cases_Cumulative) - 
         subset(cumcases.state, Date == earliestdate, select = Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
current_ratio[current_ratio<0] <- NA

state_current_ratio_dat <- 
  data.frame(current_ratio=current_ratio,
             current_ratio_cat = cut(current_ratio,
                                     breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

export

write.csv(state_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/state_case_ratio.csv', row.names = F)
write.csv(tsa_current_ratio_dat,'statistical-output/standard-stats/case-ratios/tsa_case_ratio.csv', row.names = F)
write.csv(county_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/county_case_ratio.csv', row.names = F)
write.csv(metro_current_ratio_dat,'statistical-output/standard-stats/case-ratios/metro_case_ratio.csv', row.names = F)
write.csv(phr_current_ratio_dat,'statistical-output/standard-stats/case-ratios/phr_case_ratio.csv', row.names = F)

write.csv(state_current_ratio_dat, 'tableau/case-ratios/state_case_ratio.csv', row.names = F)
write.csv(tsa_current_ratio_dat,'tableau/case-ratios/tsa_case_ratio.csv', row.names = F)
write.csv(county_current_ratio_dat, 'tableau/case-ratios/county_case_ratio.csv', row.names = F)
write.csv(metro_current_ratio_dat,'tableau/case-ratios/metro_case_ratio.csv', row.names = F)
write.csv(phr_current_ratio_dat,'tableau/case-ratios/phr_case_ratio.csv', row.names = F)

grouping

colnames(county_current_ratio_dat)[1] <- 'Level'
colnames(metro_current_ratio_dat)[1] <- 'Level'
colnames(tsa_current_ratio_dat)[1] <- 'Level'
colnames(phr_current_ratio_dat)[1] <- 'Level'
state_current_ratio_dat$Level <- 'Texas'


county_current_ratio_dat$Level_Type = 'County'
metro_current_ratio_dat$Level_Type = 'Metro'
tsa_current_ratio_dat$Level_Type = 'TSA'
phr_current_ratio_dat$Level_Type = 'PHR'
state_current_ratio_dat$Level_Type = 'State'

combined_current_ratio_dat <- rbind(county_current_ratio_dat, metro_current_ratio_dat,
                                   tsa_current_ratio_dat, state_current_ratio_dat, phr_current_ratio_dat)

write.csv(combined_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/stacked_case_ratio.csv',
          row.names = FALSE)

# write.csv(combined_current_ratio_dat, 'tableau/stacked_case_ratio.csv',
#           row.names = FALSE)

% CHANGE

create.case.test <-   function(level, dat, region){
  # creates the % difference in cases and tests and smooth line with CIs 
  # level: either "TSA", "county", or "metro". Note that "county" won't work for many counties unless have enough cases.
  # dat: dataset (e.g. "county", "metro", "tsa")
  # region: the region within the dataset (county, metro region, or tsa)
  
  if(level == "TSA"){
    dat <- subset(dat, TSA==region)
  }
   if(level == "PHR"){
    dat <- subset(dat, PHR==region)
  }
  if(level == "county"){
    dat <- subset(dat, County == region)
  }
  if(level == "metro"){
    dat <- subset(dat, Metro_Area==region)
  }

  # restrict data to first test date (for % test increase)
  first.test.date <- ymd("2020-04-22")
  dat$Date <- ymd(dat$Date)
  dat$Date2 <- as.numeric(ymd(dat$Date))
  
  # get 14 day rolling avg for instances where initial cases or tests are 0
  dat$cases_avg14 <- rollmean(dat$Cases_Daily, k=14, fill = NA, align = 'right', na.rm = TRUE)
  dat$tests_avg14 <- rollmean(dat$Tests_Daily, k=14, fill = NA, align = 'right', na.rm = TRUE)
  
  # restrict new df to first.test.date forward  
  slopedata.tests <- data.frame(subset(dat, select=c(Date, Date2, Cases_Daily,
                                                     Tests_Daily, cases_avg14, tests_avg14))) %>%
    subset(ymd(Date) >= ymd(first.test.date))
  
  tests_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=Tests_Daily))
  cases_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=Cases_Daily))
  
  # numeric date for gam                         
  tmpdata <- data.frame(Date2=slopedata.tests$Date2)

  if (cases_start != 0 & tests_start != 0) {
    
    # Calc splines if cases and tests will produce non Inf values
    slopedata.tests$newtestsY <- 100*(as.vector(slopedata.tests$Tests_Daily / tests_start) - 1)
    slopedata.tests$newcasesY <- 100*(as.vector(slopedata.tests$Cases_Daily / cases_start) - 1)
    
    # fit and predict w/ gam for spline vals
    cases.gam <- gam(newcasesY~s(Date2,bs="cs", k=5), data=slopedata.tests)
    tests.gam <- gam(newtestsY~s(Date2,bs="cs",k=5), data=slopedata.tests)
    
    cases.fit <- predict(cases.gam, tmpdata, se.fit=TRUE)
    tests.fit <- predict(tests.gam, tmpdata, se.fit=TRUE)
    
    # Use 95% CI for splines
    tmp.df <- data.frame(Date = slopedata.tests$Date,
                         Data_Type = 'Raw',
                         cases_percentdiff = slopedata.tests$newcasesY,
                         tests_percentdiff = slopedata.tests$newtestsY,
                         cases_percentdiff_spline = cases.fit$fit, 
                         cases_percentdiff_spline_lower = cases.fit$fit-1.96*cases.fit$se, 
                         cases_percentdiff_spline_upper = cases.fit$fit+1.96*cases.fit$se,
                         tests_percentdiff_spline = tests.fit$fit,
                         tests_percentdiff_spline_lower = tests.fit$fit-1.96*tests.fit$se,
                         tests_percentdiff_spline_upper = tests.fit$fit+1.96*tests.fit$se)    
  } else {
    
    # only calc % increase for regions w/ sparse cases
    cases_avg_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=cases_avg14))
    tests_avg_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=tests_avg14))

    slopedata.tests$newcasesY <- 100*(as.vector(slopedata.tests$cases_avg14 / cases_avg_start) - 1)
    slopedata.tests$newtestsY <- 100*(as.vector(slopedata.tests$tests_avg14 / tests_avg_start) - 1)
    
    tmp.df <- data.frame(Date = slopedata.tests$Date,
                         Data_Type = '14-Day Average',
                         cases_percentdiff = slopedata.tests$newcasesY,
                         tests_percentdiff = slopedata.tests$newtestsY,
                         cases_percentdiff_spline = rep(NA, times = nrow(slopedata.tests)),
                         cases_percentdiff_spline_lower = rep(NA, times = nrow(slopedata.tests)),
                         cases_percentdiff_spline_upper = rep(NA, times = nrow(slopedata.tests)),
                         tests_percentdiff_spline = rep(NA, times = nrow(slopedata.tests)),
                         tests_percentdiff_spline_lower = rep(NA, times = nrow(slopedata.tests)),
                         tests_percentdiff_spline_upper = rep(NA, times = nrow(slopedata.tests)))  
    }
  
  return(tmp.df)  
}

State

state <- readxl::read_xlsx('combined-datasets/state.xlsx', sheet=1)
state.case.test.df <- create.case.test(level="State", state, NA)
write.csv(state.case.test.df, 'statistical-output/standard-stats/pct-change/state_pct_change.csv',
          row.names = FALSE)

TSA

tsa <- read.csv('combined-datasets/tsa.csv')

# Create tsa-level data (can optimize this code in the future, but it runs pretty quickly already)
# TODO: convert to gapply -> rbindlist
tsalist <- unique(tsa$TSA)
tmp <- create.case.test(level="TSA", tsa,tsalist[1])
tsa.case.test.df <- data.frame(TSA=rep(tsalist[1], nrow(tmp)), 
                               TSA_Name=rep(unique(tsa$TSA_Name[1])),
                               tmp)

for(i in c(2:length(tsalist))){
  tmp <- create.case.test(level="TSA", tsa,tsalist[i])
  tsa.case.test.df <- rbind(tsa.case.test.df, data.frame(TSA=rep(tsalist[i], nrow(tmp)), 
                                                         TSA_Name=rep(unique(tsa$TSA_Name[i])),
                                                         tmp))
}


## combine tsa name
tsa.case.test.df$TSA = paste0(tsa.case.test.df$TSA, ' - ', tsa.case.test.df$TSA_Name)
tsa.case.test.df$TSA_Name = NULL

write.csv(tsa.case.test.df, 'statistical-output/standard-stats/pct-change/tsa_pct_change.csv',
          row.names = FALSE)

PHR

phr <- read.csv('combined-datasets/phr.csv')

# Create phr-level data (can optimize this code in the future, but it runs pretty quickly already)
# TODO: convert to gapply -> rbindlist
phrlist <- unique(phr$PHR)
tmp <- create.case.test(level="PHR", phr,phrlist[1])
phr.case.test.df <- data.frame(PHR=rep(phrlist[1], nrow(tmp)), 
                               PHR_Name=rep(unique(phr$PHR_Name[1])),
                               tmp)

for(i in c(2:length(phrlist))){
  tmp <- create.case.test(level="PHR", phr,phrlist[i])
  phr.case.test.df <- rbind(phr.case.test.df, data.frame(PHR=rep(phrlist[i], nrow(tmp)), 
                                                         PHR_Name=rep(unique(phr$PHR_Name[i])),
                                                         tmp))
}


## combine phr name
phr.case.test.df$PHR = paste0(phr.case.test.df$PHR, ' - ', phr.case.test.df$PHR_Name)
phr.case.test.df$PHR_Name = NULL

write.csv(phr.case.test.df, 'statistical-output/standard-stats/pct-change/phr_pct_change.csv',
          row.names = FALSE)

Metro

Metro <- read.csv('combined-datasets/metro.csv')

metrolist <- unique(Metro$Metro_Area)

tmp <- create.case.test(level="metro", Metro,metrolist[1])
metro.case.test.df <- data.frame(Metro_Area=rep(metrolist[1], nrow(tmp)), tmp)

tmp <- create.case.test(level="metro", Metro,metrolist[2])
metro.case.test.df <- rbind(metro.case.test.df, data.frame(Metro_Area=rep(metrolist[2], nrow(tmp)), tmp))

write.csv(metro.case.test.df, 'statistical-output/standard-stats/pct-change/metro_pct_change.csv',
          row.names = FALSE)

County

county <- read.csv('combined-datasets/county.csv')
countylist <- unique(county$County)

tmp <- create.case.test(level="county", county,countylist[1])
county.case.test.df <- data.frame(County=rep(countylist[1], nrow(tmp)), tmp)

for(i in 2:length(countylist)) {
  tmp <- create.case.test(level="county", county,countylist[i])
  county.case.test.df <- rbind(county.case.test.df, data.frame(County=rep(countylist[i], nrow(tmp)), tmp))
}

# assess missing vals
missing_vals <- sum(is.na(county.case.test.df$cases_percentdiff)) + sum(is.na(county.case.test.df$tests_percentdiff))
recorded_vals <- sum(!is.na(county.case.test.df$cases_percentdiff)) + sum(!is.na(county.case.test.df$tests_percentdiff))

# 11.88% w/ 2020-05-01
# 12.71% w/ 2020-04-22
# 12% w/ 2020-05-15
missing_vals / (missing_vals + recorded_vals)
## [1] 0.1016765
write.csv(county.case.test.df, 'statistical-output/standard-stats/pct-change/county_pct_change.csv',
          row.names = FALSE)

grouping

colnames(county.case.test.df)[1] <- 'Level'
colnames(metro.case.test.df)[1] <- 'Level'
colnames(tsa.case.test.df)[1] <- 'Level'
colnames(phr.case.test.df)[1] <- 'Level'
state.case.test.df$Level <- 'Texas'

county.case.test.df$Level_Type = 'County'
metro.case.test.df$Level_Type = 'Metro'
tsa.case.test.df$Level_Type = 'TSA'
phr.case.test.df$Level_Type = 'PHR'
state.case.test.df$Level_Type = 'State'

combined.case.test.df <- rbind(county.case.test.df, metro.case.test.df,
                               tsa.case.test.df, state.case.test.df, phr.case.test.df)
write.csv(combined.case.test.df, 'statistical-output/standard-stats/pct-change/stacked_pct_change.csv',
          row.names = FALSE)

write.csv(combined.case.test.df, 'tableau/stacked_pct_change.csv',
          row.names = FALSE)

STACK COMBINATIONS

combined_current_ratio_dat$Date = max(combined.case.test.df$Date)
colnames(combined.arima.df.1)[5:6] = c('TS_CI_Lower', 'TS_CI_Upper')
colnames(combined.rt.df)[4:5] = c('RT_CI_Lower', 'RT_CI_Upper')

stacked_all = Reduce(function(x, y) merge(x, y, by = c('Level_Type', 'Level', 'Date'), all=TRUE),
       list(combined.case.test.df, combined_current_ratio_dat, combined.arima.df.1, combined.rt.df))

write.csv(stacked_all, 'tableau/stacked_all.csv', row.names = FALSE)

rolling average

Omitted fow now, Dr. Yaseen is processing this in Tableau